discover and transform plastic waste makers index data

purpose of notebook

(-) test and play with advanced numerical EDA methods

todos: (-) …

information

name: makeovermonday_2021w22 link: https://data.world/makeovermonday/2021w22 title: 2021/W22: The Plastic Waste Makers Index Data Source: Minderoo from 2019

insights

1. correlation - most of the columns are highly correlated, that was to be expected, since most variables are depend on each other, e.g., rigid + flexible = total c production, total -> -rank rigid overall has a less strong correlation with the other variables, which might hint to there being a a different sub-population based on rigid production, flexible has a stronger correlation with total as rigid and total, since flexible is a far bigger contribution to total flexible has a stronger correlation with the overall production, than rigid, this is interesting and also might hint to rigid producers have a different market strategy than flexible producers
2. a Gini coefficient of 0.56 is quite high, which means that only a top few producers are responsible for a large amount of the total SUP waste contribution, in other words 10% are responsible for 44% of the waste
3. I need more knowledge how to work with and interpret SOM and PCA, maybe also not enough observations in data set,
   but the clustering from SOM, when only independent variables are used, actually show the two assumed subpopulations based on rigid vs flexible production
4. parallel coordinate plot is not very insightful, since there is so much data close together, we cannot see any grouping without prior knowledge, only hint is when you zoom in at the bulk of the data in the lower third (y direction), there is a negative correlation
5. remind that the difference with the partition by k-means is that for hierarchical clustering, the number of classes is not specified in advance
   it seems that most clusters are confused by some specific values (81, 63, 72, 36?), are they some sort of outlier, or what is special about them?
6. As a reminder, this method aims at partitioning n observations into k clusters in which each observation belongs to the cluster with the closest average, serving as a prototype of the cluster works not quite as well as SOM, but very close, also recommend cluster is three, but closely followed by 2, since 3 does not reveal any significant connection right now, I wonder what that shows

references

1. https://iamciera.github.io/SOMexample/html/SOM_RNAseq_tutorial_part2a_SOM.html
2. https://statsandr.com/blog/clustering-analysis-k-means-and-hierarchical-clustering-by-hand-and-in-r/#optimal-number-of-clusters
3. https://lukedaniels1.github.io/Bio381_2018/Daniels_Cluster_Analysis_Lecture.html
4. https://towardsdatascience.com/10-tips-for-choosing-the-optimal-number-of-clusters-277e93d72d92

load packages & global options

overview

library(tidyverse) # tidy data frame
library(ggthemes) # for extra plot themes
library(plotly) # make ggplots interactive

library(factoextra) # provides some easy-to-use functions to extract and visualize the output of multivariate data analyses

# individual libraries are in the according cell

knitr::opts_chunk$set(
  # fig.width = 15, fig.height = 9, 
  warning = FALSE
)

# plotly: ,width = 900, height = 550

head(plastic)

correlation
most of the columns are highly correlated, that was to be expected, since most variables are depend on each other, e.g., rigid + flexible = total c production, total -> -rank rigid overall has a less strong correlation with the other variables, which might hint to there being a a different sub-population based on rigid production, flexible has a stronger correlation with total as rigid and total, since flexible is a far bigger contribution to total flexible has a stronger correlation with the overall production, than rigid, this is interesting and also might hint to rigid producers have a different market strategy than flexible producers
%>% select(-rank, -total, -assets) can show a more clear picture by removing dependent variables
r summary(plastic)
polymer_producer rank no_of_assets production_of_in_scope_polymers flexible_format_contribution_to_sup_waste rigid_format_contribution_to_sup_waste Length:100 Min. : 1.00 Min. : 0.00 Min. : 0.200 Min. :0.000 Min. :0.000 Class :character 1st Qu.: 25.75 1st Qu.: 3.00 1st Qu.: 0.500 1st Qu.:0.100 1st Qu.:0.100 Mode :character Median : 50.50 Median : 6.00 Median : 0.900 Median :0.200 Median :0.200 Mean : 50.50 Mean :11.56 Mean : 1.805 Mean :0.538 Mean :0.416 3rd Qu.: 75.25 3rd Qu.:12.25 3rd Qu.: 1.700 3rd Qu.:0.500 3rd Qu.:0.500 Max. :100.00 Max. :82.00 Max. :11.600 Max. :4.700 Max. :4.500 total_contribution_to_sup_waste total_waste_div_production Min. :0.200 Min. :0.3000 1st Qu.:0.300 1st Qu.:0.4300 Median :0.450 Median :0.5000 Mean :0.950 Mean :0.5834 3rd Qu.:0.925 3rd Qu.:0.6900 Max. :5.900 Max. :1.0000

Lorenz curve & Gini coefficient

a Gini coefficient of 0.56 is quite high, which means that only a top few producers are responsible for a large amount of the total SUP waste contribution, in other words 10% are responsible for 44% of the waste

name = c('')
df <- plastic %>% select(-polymer_producer, - total_waste_div_production) %>% mutate(rank = -rank) %>% # change sign of rank to make it increase with the dependent variables
  rename( assets = no_of_assets, 
          product = production_of_in_scope_polymers, 
          flexible = flexible_format_contribution_to_sup_waste, 
          rigid = rigid_format_contribution_to_sup_waste, 
          total = total_contribution_to_sup_waste) 


library(corrplot) # correlation plots
# https://cran.r-project.org/web/packages/corrplot/vignettes/corrplot-intro.html

cor <- cor(df)
cor_mtest <- cor.mtest(df, conf.level = 0.99) # combining correlogram with significance test
corrplot(cor, method = "number", order = 'hclust', addrect = 3, p.mat = cor_mtest$p, insig = "pch") # order = AOE, FPC, hclust + addrect

corrplot(cor, p.mat = cor_mtest$p, low = cor_mtest$lowCI, upp = cor_mtest$uppCI, order = 'hclust', sig.level = 0.01, tl.pos = 'd', addrect = 3, rect.col = 'navy', plotC = 'rect', cl.pos = 'n', insig = "p-value")

name = c('total_contribution_to_sup_waste', 'Producers')
df <- plastic %>% rename(value = total_contribution_to_sup_waste) %>% select(value) 

library(gglorenz) #transformations for plotting Lorenz curve, https://github.com/jjchern/gglorenz

lorenzcurve <- df %>% 
  ggplot(aes(value)) +
    stat_lorenz(desc = FALSE) +
    coord_fixed() +
    geom_abline(linetype = 'dashed') +
    theme_minimal() +
    hrbrthemes::scale_x_percent() +
    hrbrthemes::scale_y_percent() +
    # hrbrthemes::theme_ipsum_rc() +
    annotate_ineq(df$value) +
    ggtitle(paste("Lorenz curve for", name[1], sep=" ")) 
lorenzcurve <- ggplotly(lorenzcurve) %>% layout(yaxis = list(title = paste("cumulative percentage of", name[1], sep=" ")), xaxis = list(title = paste("cumulative percentage of", name[2], sep=" "))) 

lorenzcurve

principal component analysis colored by self organizing map cluster

I need more knowledge how to work with and interpret SOM and PCA, maybe also not enough observations in data set https://iamciera.github.io/SOMexample/html/SOM_RNAseq_tutorial_part2a_SOM.html

but the clustering from SOM, when only independent variables are used, actually show the two assumed subpopulations based on rigid vs flexible produciton

name = c('flexible_format_contribution_to_sup_waste', 'rigid_format_contribution_to_sup_waste', 'Producers', 'flexible', 'rigid')
df <- plastic %>% rename(flexible = flexible_format_contribution_to_sup_waste, rigid = rigid_format_contribution_to_sup_waste) %>% select(flexible, rigid) %>% pivot_longer(cols = c(flexible,rigid))

library(gglorenz) #transformations for plotting Lorenz curve, https://github.com/jjchern/gglorenz

# get ggplot standard colors for grouping, which are equally spaced hues around the color wheel, starting from 15
gg_color_hue <- function(n) {
  hues = seq(15, 375, length = n + 1)
  hcl(h = hues, l = 65, c = 100)[1:n]
}

lorenzcurve <- df %>% 
  ggplot(aes(x = value, color = name)) +
    stat_lorenz(desc = FALSE) +
    coord_fixed() +
    geom_abline(linetype = 'dashed') +
    theme_minimal() +
    hrbrthemes::scale_x_percent() +
    hrbrthemes::scale_y_percent() +
    # hrbrthemes::theme_ipsum_rc() +
    annotate_ineq(filter(df, name == name[4])$value, y = 0.95, colour = gg_color_hue(2)[1]) +
    annotate_ineq(filter(df, name == name[5])$value, y = 0.90, colour = gg_color_hue(2)[2]) +
    ggtitle(paste("compare Lorenz curve of", name[1], "and", name[2], sep=" ")) 
lorenzcurve <- ggplotly(lorenzcurve) %>% layout(yaxis = list(title = paste("cumulative percentage of<br>", name[1], "<br>", name[2], sep="")), xaxis = list(title = paste("cumulative percentage of", name[2], sep=" "))) 

lorenzcurve

name = c('polymer_producer')
df <- plastic %>% select(- total_waste_div_production, -rank, -no_of_assets, -total_contribution_to_sup_waste) %>% # removed variables which are depended on each other
  rename( product = production_of_in_scope_polymers, 
          flexible = flexible_format_contribution_to_sup_waste, 
          rigid = rigid_format_contribution_to_sup_waste)


library(kohonen) # functions to train self-organising maps (SOMs)

# setup for pca
scale_data <- as.matrix(t(scale(t(df[, !names(df) %in% name])))) # We need to normalize the data based on scale function because the variables are different scales; Normalization means subtracting mean from each observation and dividing with standard deviation. Check all the variables mean values are zero now
head(scale_data)

       product   flexible      rigid
[1,] 1.0838622 -0.1970659 -0.8867964
[2,] 1.0301899 -0.0633963 -0.9667936
[3,] 1.1172018 -0.3058262 -0.8113756
[4,] 0.6859477 -1.1474034  0.4614557
[5,] 1.1208971 -0.3202563 -0.8006408
[6,] 1.0995250 -0.2443389 -0.8551861

# principle component analysis
pca <- prcomp(scale_data, scale=TRUE)
summary(pca)

Importance of components:
                          PC1     PC2       PC3
Standard deviation     1.6591 0.49736 4.375e-16
Proportion of Variance 0.9175 0.08246 0.000e+00
Cumulative Proportion  0.9175 1.00000 1.000e+00

# visualize pcs results
# Contributions of variables to PC1
fviz_contrib(pca, choice = "var", axes = 1, top = 10)

# Contributions of variables to PC2
fviz_contrib(pca, choice = "var", axes = 2, top = 10)

# Control variable colors using their contributions to the principle axis
fviz_pca_var(pca, col.var="contrib",
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel = TRUE # Avoid text overlapping
             ) + theme_minimal() + ggtitle("Variables - PCA")

# add back to original so everything is together
pca_scores <- data.frame(pca$x)
data_val <- cbind(df, pca_scores)

pca_plot <- ggplot(data_val, aes(x = PC1, y = PC2)) +
    geom_rug(alpha = 0.5) + # two 1d marginal distributions, display individual cases so are best used with smaller datasets
    geom_density_2d(alpha = 0.2, bins = 4) +# 2D kernel density estimation using MASS::kde2d() and display the results with contours
    geom_point(alpha = 0.75) + # point geom is used to create scatterplots
    theme_minimal()
pca_plot <- ggplotly(pca_plot) %>% layout()

pca_plot

# clustering is performed using the som() function on the scaled gene expression values.
set.seed(3)

# define a grid for the SOM and train
grid_size <- ncol(scale_data)
som_grid <- somgrid(xdim = grid_size, ydim = grid_size, topo = 'hexagonal')
som_model <- som(scale_data, grid = som_grid)
summary(som_model)

SOM of size 3x3 with a hexagonal topology and a bubble neighbourhood function.
The number of data layers is 1.
Distance measure(s) used: sumofsquares.
Training data included: 100 objects.
Mean distance to the closest unit in the map: 0.007.

# generate som plots after training
plot(som_model, type = 'mapping')

plot(som_model, type = 'codes')

# plot(som_model, type = 'counts')
# plot(som_model, type = 'dist.neighbours')
# plot(som_model, type = 'quality')
# plot(som_model, type = 'changes')

# further split the clusters into a smaller set of clusters using hierarchical clustering.
som_cluster <- cutree(hclust(dist(som_model$codes[[1]])), 2) # use hierarchical clustering to cluster the codebook vectors

plot(som_model, type="mapping", bgcol = som_cluster, main = "Clusters")
add.cluster.boundaries(som_model, som_cluster)

# attach the hierchal cluster to the larger dataset data_val.
gridSquare <- grid_size * grid_size
som_clusterKey <- data.frame(som_cluster)
som_clusterKey$unit_classif <- c(1:gridSquare)
data_val <- cbind(data_val,som_model$unit.classif,som_model$distances) %>% rename(unit_classif = 'som_model$unit.classif', distances = 'som_model$distances')
data_val <- merge(data_val, som_clusterKey, by.x = "unit_classif" )
head(data_val)

# plot pca with colored clusters
pcasom_plot <- ggplot(data_val, aes(x = PC1, y = PC2, color = factor(som_cluster), text = polymer_producer)) +
    geom_rug(alpha = 0.5) + # two 1d marginal distributions, display individual cases so are best used with smaller datasets
    geom_point(alpha = 0.75) + # point geom is used to create scatterplots
    theme_minimal()
pcasom_plot <- ggplotly(pcasom_plot) %>% layout()

pcasom_plot

parallel coordinate plot

parallel coordinate plot is not very insightful, since there is so much data close together, we cannot see any grouping without prior knowledge, only hint is when you zoom in at the bulk of the data in the lower third (y direction), there is a negative correlation

# two variables, continuous x, continuous y, show trend and distribution
name = c('production_of_in_scope_polymers', 'total_contribution_to_sup_waste')
df <- merge(plastic, data_val, by.x = 'polymer_producer')
df <- df %>% rename(x = production_of_in_scope_polymers, y = total_contribution_to_sup_waste, cluster = som_cluster, text = polymer_producer) %>% select(x, y, cluster, text) 

# https://ggplot2.tidyverse.org/reference/geom_smooth.html
point_plot <- df %>%
  ggplot(aes(x = x, y = y, color = factor(cluster))) +
    # geom_jitter(alpha = 0.5, size = 1) +
    geom_rug(alpha = 0.5) + # two 1d marginal distributions, display individual cases so are best used with smaller datasets
    geom_density_2d(alpha = 0.2, bins = 4) +# 2D kernel density estimation using MASS::kde2d() and display the results with contours
    geom_smooth(fill = "grey90") + # aids the eye in seeing patterns in the presence of overplotting
    geom_point(aes(text = text), alpha = 0.75) + # point geom is used to create scatterplots
    theme_minimal() +
    ggtitle(paste("trend of", name[2], "over", name[1], sep=" ")) 
point_plot <- ggplotly(point_plot) %>% layout(xaxis = list(showticklabels = FALSE))

`geom_smooth()` using method = 'loess' and formula 'y ~ x'

x_density_plot <- df %>%
  ggplot(aes(x = x, color = factor(cluster))) +
    stat_density(geom="line") + # draws kernel density estimate, which is a smoothed version of the histogram
    # geom_histogram(binwidth = 1) +
    theme_minimal() 
x_density_plot <- ggplotly(x_density_plot) %>% layout(yaxis = list(showticklabels = FALSE, showgrid = FALSE), xaxis = list(showticklabels = FALSE, showgrid = FALSE))

y_density_plot <- df %>%
  ggplot(aes(x = y, color = factor(cluster))) +
    stat_density(geom="line") + # draws kernel density estimate, which is a smoothed version of the histogram
    # geom_histogram(binwidth = 1) +
    coord_flip() +
    theme_minimal() 
y_density_plot <- ggplotly(y_density_plot) %>% layout(yaxis = list(showticklabels = FALSE, showgrid = FALSE), xaxis = list(showticklabels = FALSE, showgrid = FALSE))

# https://ggplot2.tidyverse.org/reference/geom_quantile.html
qualtile_plot <- df %>%
  ggplot(aes(x = x, y = y, color = factor(cluster))) +
    geom_quantile(alpha = 0.8) + # fits a quantile regression to the data and draws the fitted quantiles with lines
    theme_minimal() 
qualtile_plot <- ggplotly(qualtile_plot) %>% layout(yaxis = list(showticklabels = FALSE, showgrid = FALSE))

Smoothing formula not specified. Using: y ~ x
Smoothing formula not specified. Using: y ~ x

# merge figures into one plot, via subplots, https://plotly-r.com/arranging-views.html
sub1 <- subplot(x_density_plot, plotly_empty(), point_plot, y_density_plot, nrows = 2, margin = 0, heights = c(0.15, 0.85), widths = c(0.9, 0.1), shareX = TRUE, shareY = TRUE, titleX = FALSE, titleY = FALSE) %>% layout()

No trace type specified:
  Based on info supplied, a 'scatter' trace seems appropriate.
  Read more about this trace type -> https://plotly.com/r/reference/#scatter
No scatter mode specifed:
  Setting the mode to markers
  Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode

sub2 <- subplot(qualtile_plot, plotly_empty(), margin = 0, widths = c(0.9, 0.10), titleX = FALSE, titleY = FALSE) %>% layout()

No trace type specified:
  Based on info supplied, a 'scatter' trace seems appropriate.
  Read more about this trace type -> https://plotly.com/r/reference/#scatter
No scatter mode specifed:
  Setting the mode to markers
  Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode

fig <- subplot(sub1, sub2, nrows = 2, margin = 0, heights = c(0.8, 0.2), shareX = TRUE) %>% layout(xaxis = list(title = name[1]), yaxis = list(title = name[2]))
  
fig

k-means clustering

https://statsandr.com/blog/clustering-analysis-k-means-and-hierarchical-clustering-by-hand-and-in-r/#optimal-number-of-clusters As a reminder, this method aims at partitioning n observations into k clusters in which each observation belongs to the cluster with the closest average, serving as a prototype of the cluster works not quite as well as SOM, but very close, also recommend cluster is three, but closely followed by 2, since 3 does not reveal any significant connection right now, I wonder what that shows

name = c('plastic producers clustered by focus')
df <- merge(plastic, select(data_val, som_cluster, polymer_producer), by.x = 'polymer_producer') %>%
  mutate(som_cluster = as.character(som_cluster)) %>%
  rename(product = production_of_in_scope_polymers, 
          flexible = flexible_format_contribution_to_sup_waste, 
          rigid = rigid_format_contribution_to_sup_waste,
          assets =  no_of_assets,
          cluster = som_cluster) %>% 
  select(assets, product, flexible, rigid, cluster) 

library(GGally) # extends ggplot2 by adding several functions to reduce the complexity of combining geoms with transformed data

# https://www.r-graph-gallery.com/parallel-plot-ggally.html#custom
parcoord_plot <- ggparcoord(df,
           columns = 1:4, groupColumn = 5,
           scale='center', # scaling: standardize and center variables
           showPoints = TRUE,
           title = name,
           alphaLines = 0.3) +
  theme_minimal() 
parcoord_plot <- ggplotly(parcoord_plot) %>% layout(autosize=T)

parcoord_plot

name = c('polymer_producer')
df <- plastic %>% select(- total_waste_div_production, -rank, -no_of_assets, -total_contribution_to_sup_waste) %>% # removed variables which are depended on each other
  rename( product = production_of_in_scope_polymers, 
          flexible = flexible_format_contribution_to_sup_waste, 
          rigid = rigid_format_contribution_to_sup_waste)

library(NbClust) # determining the optimal number of clusters in a data set
library(cluster) # computes agglomerative hierarchical clustering of the dataset

scale_data <- as.matrix(t(scale(t(df[, !names(df) %in% name])))) # We need to normalize the data based on scale function because the variables are different scales; Normalization means subtracting mean from each observation and dividing with standard deviation. Check all the variables mean values are zero now
head(scale_data)

       product   flexible      rigid
[1,] 1.0838622 -0.1970659 -0.8867964
[2,] 1.0301899 -0.0633963 -0.9667936
[3,] 1.1172018 -0.3058262 -0.8113756
[4,] 0.6859477 -1.1474034  0.4614557
[5,] 1.1208971 -0.3202563 -0.8006408
[6,] 1.0995250 -0.2443389 -0.8551861

kmeans_model <- kmeans(scale_data, centers = 2, nstart = 12) #  k-means clustering is done via the kmeans() function, with the argument centers that corresponds to the number of desired clusters
df_cluster <- tibble(df, cluster = as.factor(kmeans_model$cluster)) # store cluster in original data set as column
head(df_cluster)

# check quality of a k-means partition
quality <- kmeans_model$betweenss / kmeans_model$totss 
print(paste("quality of kmeans is BSS/TSS: ", format(round(quality,2), nsmall = 2)))

[1] "quality of kmeans is BSS/TSS:  0.88"

# find optimal number of clusters
fviz_nbclust(scale_data, kmeans, method = 'wss') + # Elbow method, needs scaled data
  geom_vline(xintercept = 2, linetype = 2) + # add line for better visualization
  labs(subtitle = "Elbow method") # add subtitle

fviz_nbclust(scale_data, kmeans, method = 'silhouette') + # Silhouette method, needs scaled data
  labs(subtitle = "Silhouette method") # add subtitle

fviz_nbclust(df[, !names(df) %in% name], kmeans, # Gap statistics, needs original data ?
             nstart = 25,
             method = 'gap_stat',
             nboot = 100) + # reduce it for lower computation time, but less precise results
  labs(subtitle = "Gap statistics method")

Clustering k = 1,2,..., K.max (= 10): .. done
Bootstrapping, b = 1,2,..., B (= 100)  [one "." per sample]:
.................................................. 50 
.................................................. 100 

nbclust_out <- NbClust(data = df[, !names(df) %in% name], # NbClust needs the original data ?
                       distance = 'euclidean',
                       min.nc = 2, # minimum number of clusters
                       max.nc = 10, # maximum number of cluster
                       method = 'complete',
                       index = 'all')

[1] "Frey index : No clustering structure in this data set"
*** : The Hubert index is a graphical method of determining the number of clusters.
                In the plot of Hubert index, we seek a significant knee that corresponds to a 
                significant increase of the value of the measure i.e the significant peak in Hubert
                index second differences plot. 
 

*** : The D index is a graphical method of determining the number of clusters. 
                In the plot of D index, we seek a significant knee (the significant peak in Dindex
                second differences plot) that corresponds to a significant increase of the value of
                the measure. 
 
******************************************************************* 
* Among all indices:                                                
* 7 proposed 2 as the best number of clusters 
* 9 proposed 3 as the best number of clusters 
* 1 proposed 5 as the best number of clusters 
* 2 proposed 6 as the best number of clusters 
* 2 proposed 9 as the best number of clusters 
* 2 proposed 10 as the best number of clusters 

                   ***** Conclusion *****                            
 
* According to the majority rule, the best number of clusters is  3 
 
 
******************************************************************* 

fviz_nbclust(nbclust_out) + theme_minimal() +
  labs(subtitle = "NbClust results")

Among all indices: 
===================
* 2 proposed  0 as the best number of clusters
* 7 proposed  2 as the best number of clusters
* 9 proposed  3 as the best number of clusters
* 1 proposed  5 as the best number of clusters
* 2 proposed  6 as the best number of clusters
* 2 proposed  9 as the best number of clusters
* 2 proposed  10 as the best number of clusters
* 1 proposed  NA's as the best number of clusters

Conclusion
=========================
* According to the majority rule, the best number of clusters is  3 .

# check quality of clustering
# if a large majority of the silhouette coefficients are positive, it indicates that the observations are placed in the correct group
sil <- silhouette(kmeans_model$cluster, dist(scale_data)) 
fviz_silhouette(sil)

fviz_cluster(kmeans_model, df[, !names(df) %in% name], ellipse.type = 'norm') + theme_minimal()

name = c('polymer_producer')
df <- plastic %>% select(- total_waste_div_production, -rank, -no_of_assets, -total_contribution_to_sup_waste) %>% # removed variables which are depended on each other
  rename( product = production_of_in_scope_polymers, 
          flexible = flexible_format_contribution_to_sup_waste, 
          rigid = rigid_format_contribution_to_sup_waste)

scale_data <- as.matrix(t(scale(t(df[, !names(df) %in% name])))) # We need to normalize the data based on scale function because the variables are different scales; Normalization means subtracting mean from each observation and dividing with standard deviation. Check all the variables mean values are zero now

kmean_calc <- function(df, ...){
  kmeans(df, scaled = ..., nstart = 30)
}

km2 <- kmean_calc(scale_data, 2)
km3 <- kmean_calc(scale_data, 3)
km4 <- kmeans(scale_data, 4)
km5 <- kmeans(scale_data, 5)
km6 <- kmeans(scale_data, 6)
km7 <- kmeans(scale_data, 7)

p1 <- fviz_cluster(km2, data = scale_data, ellipse.type = "convex") + theme_minimal() + ggtitle("k = 2") 
p2 <- fviz_cluster(km3, data = scale_data, ellipse.type = "convex") + theme_minimal() + ggtitle("k = 3")
p3 <- fviz_cluster(km4, data = scale_data, ellipse.type = "convex") + theme_minimal() + ggtitle("k = 4")
p4 <- fviz_cluster(km5, data = scale_data, ellipse.type = "convex") + theme_minimal() + ggtitle("k = 5")
p5 <- fviz_cluster(km6, data = scale_data, ellipse.type = "convex") + theme_minimal() + ggtitle("k = 6")
p6 <- fviz_cluster(km7, data = scale_data, ellipse.type = "convex") + theme_minimal() + ggtitle("k = 7")

fig <- subplot(p1, p2, p3 , p4, p5, p6, nrows = 2, shareX = TRUE, shareY = TRUE) %>% layout(autosize=T) # TOOD: make all plots linked
fig

hierarchical clustering

https://statsandr.com/blog/clustering-analysis-k-means-and-hierarchical-clustering-by-hand-and-in-r/#optimal-number-of-clusters remind that the difference with the partition by k-means is that for hierarchical clustering, the number of classes is not specified in advance

it seems that most clusters are confused by some specific values (81, 63, 72, 36?), are they some sort of outlier, or what is special about them?

# two variables, continuous x, continuous y, show trend and distribution
name = c('production_of_in_scope_polymers', 'total_contribution_to_sup_waste')
df <- tibble(plastic, cluster = as.factor(kmeans_model$cluster))
df <- df %>% rename(x = production_of_in_scope_polymers, y = total_contribution_to_sup_waste, cluster = cluster, text = polymer_producer) %>% select(x, y, cluster, text) 

# https://ggplot2.tidyverse.org/reference/geom_smooth.html
point_plot <- df %>%
  ggplot(aes(x = x, y = y, color = factor(cluster))) +
    # geom_jitter(alpha = 0.5, size = 1) +
    geom_rug(alpha = 0.5) + # two 1d marginal distributions, display individual cases so are best used with smaller datasets
    geom_density_2d(alpha = 0.2, bins = 4) +# 2D kernel density estimation using MASS::kde2d() and display the results with contours
    geom_smooth(fill = "grey90") + # aids the eye in seeing patterns in the presence of overplotting
    geom_point(aes(text = text), alpha = 0.75) + # point geom is used to create scatterplots
    theme_minimal() +
    ggtitle(paste("trend of", name[2], "over", name[1], sep=" ")) 
point_plot <- ggplotly(point_plot) %>% layout(xaxis = list(showticklabels = FALSE))

`geom_smooth()` using method = 'loess' and formula 'y ~ x'

x_density_plot <- df %>%
  ggplot(aes(x = x, color = factor(cluster))) +
    stat_density(geom="line") + # draws kernel density estimate, which is a smoothed version of the histogram
    # geom_histogram(binwidth = 1) +
    theme_minimal() 
x_density_plot <- ggplotly(x_density_plot) %>% layout(yaxis = list(showticklabels = FALSE, showgrid = FALSE), xaxis = list(showticklabels = FALSE, showgrid = FALSE))

y_density_plot <- df %>%
  ggplot(aes(x = y, color = factor(cluster))) +
    stat_density(geom="line") + # draws kernel density estimate, which is a smoothed version of the histogram
    # geom_histogram(binwidth = 1) +
    coord_flip() +
    theme_minimal() 
y_density_plot <- ggplotly(y_density_plot) %>% layout(yaxis = list(showticklabels = FALSE, showgrid = FALSE), xaxis = list(showticklabels = FALSE, showgrid = FALSE))

# https://ggplot2.tidyverse.org/reference/geom_quantile.html
qualtile_plot <- df %>%
  ggplot(aes(x = x, y = y, color = factor(cluster))) +
    geom_quantile(alpha = 0.8) + # fits a quantile regression to the data and draws the fitted quantiles with lines
    theme_minimal() 
qualtile_plot <- ggplotly(qualtile_plot) %>% layout(yaxis = list(showticklabels = FALSE, showgrid = FALSE))

Smoothing formula not specified. Using: y ~ x
Smoothing formula not specified. Using: y ~ x

# merge figures into one plot, via subplots, https://plotly-r.com/arranging-views.html
sub1 <- subplot(x_density_plot, plotly_empty(), point_plot, y_density_plot, nrows = 2, margin = 0, heights = c(0.15, 0.85), widths = c(0.9, 0.1), shareX = TRUE, shareY = TRUE, titleX = FALSE, titleY = FALSE) %>% layout()

No trace type specified:
  Based on info supplied, a 'scatter' trace seems appropriate.
  Read more about this trace type -> https://plotly.com/r/reference/#scatter
No scatter mode specifed:
  Setting the mode to markers
  Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode

sub2 <- subplot(qualtile_plot, plotly_empty(), margin = 0, widths = c(0.9, 0.10), titleX = FALSE, titleY = FALSE) %>% layout()

No trace type specified:
  Based on info supplied, a 'scatter' trace seems appropriate.
  Read more about this trace type -> https://plotly.com/r/reference/#scatter
No scatter mode specifed:
  Setting the mode to markers
  Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode

fig <- subplot(sub1, sub2, nrows = 2, margin = 0, heights = c(0.8, 0.2), shareX = TRUE) %>% layout(xaxis = list(title = name[1]), yaxis = list(title = name[2]))

fig

Clustree
https://towardsdatascience.com/10-tips-for-choosing-the-optimal-number-of-clusters-277e93d72d92 In this figure the size of each node corresponds to the number of samples in each cluster, and the arrows are coloured according to the number of samples each cluster receives. A separate set of arrows, the transparent ones, called the incoming node proportion, are also coloured and shows how samples from one group end up in another group — an indicator of cluster instability.
```r name = c(‘polymer_producer’) df <- plastic %>% select(- total_waste_div_production, -rank, -no_of_assets, -total_contribution_to_sup_waste) %>% # removed variables which are depended on each other rename( product = production_of_in_scope_polymers, flexible = flexible_format_contribution_to_sup_waste, rigid = rigid_format_contribution_to_sup_waste)
scale_data <- as.matrix(t(scale(t(df[, !names(df) %in% name])))) # We need to normalize the data based on scale function because the variables are different scales; Normalization means subtracting mean from each observation and dividing with standard deviation. Check all the variables mean values are zero now head(scale_data) ```
product flexible rigid [1,] 1.0838622 -0.1970659 -0.8867964 [2,] 1.0301899 -0.0633963 -0.9667936 [3,] 1.1172018 -0.3058262 -0.8113756 [4,] 0.6859477 -1.1474034 0.4614557 [5,] 1.1208971 -0.3202563 -0.8006408 [6,] 1.0995250 -0.2443389 -0.8551861
```r no_k = 2;
# Hierarchical clustering: single linkage hclust_res <- hclust(dist(scale_data), method = ‘single’) fviz_dend(hclust_res, k = no_k, rect = TRUE) ```
`guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.row names were found from a short variable and have been discarded
r # Hierarchical clustering: complete linkage hclust_res <- hclust(dist(scale_data), method = 'complete') plot(hclust_res) rect.hclust(hclust_res, k = no_k, border = 'blue')
```r
# Hierarchical clustering: average linkage hclust_res <- hclust(dist(scale_data), method = ‘average’) plot(hclust_res) rect.hclust(hclust_res, k = no_k, border = ‘blue’) ```
```r
# Hierarchical clustering: ward hclust_res <- hclust(dist(scale_data), method = ‘ward.D’) fviz_dend(hclust_res, k = no_k, rect = TRUE) ```
`guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.row names were found from a short variable and have been discarded
r # Hierarchical clustering: mcquitty hclust_res <- hclust(dist(scale_data), method = 'mcquitty') plot(hclust_res) rect.hclust(hclust_res, k = no_k, border = 'blue')
```r
# Hierarchical clustering: centroid hclust_res <- hclust(dist(scale_data), method = ‘centroid’) plot(hclust_res) rect.hclust(hclust_res, k = no_k, border = ‘blue’) ```

cValid to choose best clustering algo

https://towardsdatascience.com/10-tips-for-choosing-the-optimal-number-of-clusters-277e93d72d92 The cValid package can be used to simultaneously compare multiple clustering algorithms, to identify the best clustering approach and the optimal number of clusters. We will compare k-means, hierarchical and PAM clustering.